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ABSTRACT 



The problem of determining the acoustic pressure in a wedge shaped wave guide is 
examined. Two computer models, one using the method of images, the other using coupled 
normal modes, are compared. This comparison is over the benchmark wedge of the 
Acoustic Society of America. Three scenarios were examined: isovelocity water over a 
pressure release bottom, isovelocity water over a penetrable, lossless bottom, and isovelocity 
water over a penetrable, lossy bottom. In all cases good agreement was seen between both 
models, with some differences due to a rigid sub-bottom in the normal mode model. The 
strengths and weaknesses of each model is examined. An analytic solution in normal modes 
of the waveguide with a pressure release surface, rigid bottom, and a discontinuous speed of 
sound profile is presented. A brief history of research into the wedge problem is included in 
the introduction. 
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I. INTRODUCTION 



With the recent decline in the threat of deep water submarine forces, growing emphasis 
is being placed on operations in relatively shallow water. In contrast to open ocean models, 
shallow water models of sound propagation are influenced by interaction with the bottom. 
Variables such as bottom composition, slope angle, sediment density, etc., make acoustic 
predictions difficult at best. At present there is no closed analytic solution to the problem of 
a penetrable, sloping bottom as found close to shorelines or along the continental shelf. 
However, computer models using various approximating methods are in use. Computer 
models are limited only by the speed and memory of the computer. With faster computers 
being developed, more computationally intensive models become practical. These models 
include the parabolic equation approximation, image theory, and adiabatic normal mode 
theory. 

A. PARABOLIC EQUATION 

By replacing the Helmholz equation with a one-way parabolic equation (PE), an 
acoustic field can be generated as an initial value problem. The PE method was developed 
for the radar community and first applied to underwater acoustics by Tappert with 
restrictions on the maximum angular aperture at the source [Ref 1]. Collins finally 
eliminated this limitation of the PE model with a higher order parabolic equation [Ref 2]. PE 
methods by Jensen and Kuperman [Ref 3] took advantage of improved computing power 
and showed valid results when compared to model tank experiments [Ref 4] for upslope 
propagation. This model showed modal cutoff as the sound traveled upslope. It also showed 
the sound propagating into the bottom at the modal cutoff depths. Jensen and Tindle used 
improved numerical methods to further study this problem [Ref 5]. 
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To increase accuracy, Lee and McDaniel broke the wedge up into a number of discrete 
sections and applied the normal PE methods [Ref 6]. A more recent method of rotated 
parabolic equations by Collins has shown good results [Ref 7]. Finally, Fawcett has 
developed a three-dimensional computer model using PE methods combined with Fast 
Fourier Transform (FFT) algorithms [Ref 8]. 

B. METHOD OF IMAGES 

Similar to the image reflections in a kaleidoscope, this method combines the 
contributions of the virtual sound sources to calculate the pressure field. Macpherson and 
Dainteth proposed a phase incoherent model for upslope propagation in 1966 using the 
method of images with some success [Ref 9]. Coppens et. al. [Ref 10] later developed a 
phase coherent model to calculate the pressure amplitude along the bottom of a wedge of 
water overlying a fast fluid bottom in the upslope direction using the incident plane wave 
reflection coefficient stated by Brekhovskikh and Lysenov [Ref 1 1]. Baek then developed a 
computer model of pressure throughout the water over the fast bottom [Ref 12]. However, 
this model was limited to a fast bottom due to the critical grazing angle of the sound to the 
bottom. LeSesne used the same model as a basis for a three dimensional model, which 
produced pressure fields in the cross slope case [Ref 13], He further validated his computer 
results with experimental results from a model tank. Kaswandi then developed the model for 
a slow bottom in the downslope direction [Ref 14], Concurrently, Li introduced a model 
which applies to both fast and slow bottoms [Ref 15], Nassopoulos then took all the 
previous models underdevelopment at the Naval Postgraduate School and combined them 
in the computer program URTEXT [Ref 16] Independent of this research, Buckingham and 
Deane developed a similar model, which produced a transmission loss plot as a function of 
depth and range from the source [Ref 17]. 
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The analytic development of the method is an ongoing process. Nassopoulos 
implemented the use of image doublets [Ref 16]. Livingood then developed the image 
doublets to the three dimensional case [Ref 18]. Presently, work on Taylor expansions of 
sets of four images is being done by Joyce [Ref 19]. 

Tests done by Kim [Ref 20] and by LeSesne are cited to prove the validity of the image 
model. LeSesne compared his model with experimental conditions, while Kim compared his 
to a parabolic equation model developed by Jaeger [Ref 21] which used an implicit finite- 
difference solution. This turned out to be a favorable comparison. 

While validated in comparisons with both other methods and with experiment, the 
method of images still has some flaws. It is very computer intensive, especially as the angle 
of the bottom becomes very shallow. Improvements include truncation of images that do not 
significantly contribute to the total pressure field. 

C. NORMAL MODES 

Normal mode theory is a range independent approach. However, for the wedge 
problem a range dependent model is needed. Pierce showed that an approximation of the 
normal modes can be used by performing an adiabatic separation of depth and range in the 
wave equation [Ref 22]. Graves, Nagl, Uberall, and Zarur then applied this method to the 
wedge problem using isospeed water over a rigid bottom [Ref 23]. This method achieved 
good agreement at large angles when compared to experimental values [Ref 24], This 
approach was limited to large angles and a perfectly reflecting bottom. Buckingham 
developed a solution of the penetrable bottom. In this method he first calculates an effective 
pressure release bottom below the actual penetrable bottom. The depth of this bottom is a 
function of the density and speed of sound of the actual bottom. He then applied the sum of 
the normal modes to this new bottom [Ref 25]. 
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A problem of the adiabatic normal mode approximation is that it does not adequately 
explain the transition to evanescent modes at cutoff depths [Ref 3]. By combining PE with 
adiabatic mode theory, a critical depth function was developed which agreed with the results 
of Jensen and Kuperman [Ref 26]. 

Another method is to combine ray acoustics with local modes [Ref 27], which leads to 
the development of intrinsic wave modes [Ref 28]. Also using beam displacement with ray 
theory has shown agreement with the two-way coupled mode solution in two dimensions 
over a penetrable lossy bottom [Ref 29]. 

To take into account coupling between the modes, a version of stepwise depth 
variations was applied to normal mode theory' [Ref 30], This method w as used first by Chiu, 
Desaubies, and Miller [Ref 31], then by Chiu and Ehret into the development of a model of 
coupled normal modes [Ref 32]. The three dimensional analysis by Sagos [Ref 33] is the 
extent of published information on coupled normal modes, but improvements to the theory 
and the model are still being pursued. 
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II. DEVELOPMENT 



A. IMAGE THEORY 

In 1978, Coppens, Sanders, Iounnou, and Kawamura predicted and measured the 
pressure field and phase up slope along a wedge shaped bottom [Ref 10]. LeSesne 
implemented a model that was not dependent on up or down slope direction [Ref 13]. 
Further research cited earlier developed cross slope programs, giving three dimensional 
capability to the model. The theory behind these models is the method of images. 

Coppens et. al. provides the complete explanation of the theory but a summary is in 
order. In transmission within a waveguide there is a number of interactions of acoustic rays 
with the guide boundaries [Ref 34, p. 427]. In image theory each ray path is replaced by an 
image of the source, the distance of the source from the receiver being the total distance of 
the actual acoustic ray. For a sloping bottom the image is placed a distance from the 
reflecting boundary' equal to the distance of the source from the boundary perpendicular to 
that boundary. Figure 1 illustrates this geometry. 

The phase of the image is 180° out of phase from the source if there is an interaction 
with a pressure release surface, as in the case of the air-water interface, and in phase for a 
rigid surface— a hard bottom, for example. Each interaction of the sound with the 
boundaries implies the existence of a new image. The final pressure at the receiver can be 
calculated by taking the sum of the contributions from each image. 
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Figure 1 - Two Dimensional Wedge Geometry 

For a given bottom slope angle (B, the number of images N in each of the upper and 
lower half spaces is given by 

N='™ 

P (1) 

where p is an integer factor of 180°. Each image is numbered from the source, with the 
source being image number 1 for the upper half space. For the lower half space the first 
image below the bottom is image number 1. 

Next, the grazing angle 0„ of the «th image above or below the bottom is calculated 

from 
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e n = "P-v 



for n odd 



(2a) 



9„ = (fi - 1)P + Y for n even 

where y is the angle of the source from the surface. 

The range R n from each image to the receiver is then calculated using 

K = V r i +r 2~ 2r \ r 2 cos ( 0 « - P + 8) 

for the upper series of images and 

K = $ + r\- 2 r x r 2 cos(0„ + p - 6) 



(2b) 



(3a) 



(3b) 



for the lower series. In these equations, r\ is the range from the wedge apex to the source, ro 
is the range from the apex to the receiver, and 6 is the angle from the surface to the receiver. 

For each interaction with the surface, the reflection coefficient is -1. Each interaction 
with the bottom requires calculation of the reflection coefficient. This coefficient is a 
function of the speed of sound in each medium, the density of each medium, and the grazing 
angle of the ray equivalent on the bottom. It is sufficient to find the sine of this angle 0,„„. 
For this 



sin0„„, = 

nm 



r,sin(0„ - 2/«P) + r 2 sin[(2/n - 1$ +6] 

K 



for the upper images and 

r,sin(0 n - 2wp) + r 2 sin[(2w + l)|3 - 6] 



sin0__ - 



K 



(4a) 



(4b) 



for the lower images [Ref 16]. Figure 2 illustrates. 
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Figure 2 - Incident Angle Q nm Calculation (Upper Half Space) 

The reflection coefficient for each interaction R w „ then is given by 

P 2 _ • r\ 1 / /,2 2 • 1 / /.I -> 



R„ m = 

nm 



^ S[nQ nm-^=[^b 2 +fl2 + «) 2 +y'^(^ 2 +fl 2 -«)' 
— sin 6 „ m + - 7=(^ 2 + « 2 + «) /2 -;-^(^ 2 + « 2 



P, ^ 1 ' (5) 

where pi, p 2 are densities of the respective media,,/' W-l, and a and b are given by 



a = — I -cos 2 9, 



6 = 21 — ] a 



c 2 ) ^ 



(6a) 



(6b) 



For the above, Cj, C 2 are the speeds of sound in the respective media and a/k 2 is a 
measure of the bottom loss coefficient. [Ref 11]. We further define k n as the wave number 
in the //th medium, where 

k n = % (6c) 
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The pressure for each half space is given by summing the contribution for the upper 



and lower images respectively 

N . ... M 

p «“EF exp ("^^ (_i)W ri R ' 

n- 1 “n m ■ 1 

where ]~Jr = 1 when n = 1,2; in (7a) and 

n-l^n wO 



(7a) 



(7b) 



The number of interactions m with the bottom is related to the image number, so the limit M 
is the integer of - - * - 1 for the /zth image. For convenience we have omitted the factor 



e jml which would otherwise appear in (7a) and (7b). 

Two special cases exist in the upper half space. For the case of the direct path of the 
source to the receiver (/z=l), R is 1. In the case of a single reflection off the surface and no 
reflections with the bottom (/z=2), R is again 1. The phase inversion of the ray is computed 
by the (-1)^ term of (7a). 

Since the model is frequency independent, we will scale the ranges to a distance X: 
from the apex where the bottom depth is equal to the lowest normal mode cutoff depth when 
the bottom is fast (i.e., when the sound speed ratio is less than 1). This cutoff distance is 
calculated by 



k i x c = 



x 



\ 



0„ = cos 



-1 



\ c l) 



(8a) 



2sin0 c tan fi 

For a slow bottom an analogous scaling, convenient for the computer program is 



k x X c = 



JT 



0. = COS 



-1 



\ C 1, 



(8b) 



2 tan 6 s lan fi 

For the total pressure at the receiver, the upper and lower contributions are combined, 



thus 
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(9) 



P= P H+ P, 

B. COUPLED MODE THEORY 

Normal mode theory describes sound propagation as a collection of eigenfunctions, 
called normal modes, which are determined by the depth of the sound channel boundaries 
and the source depth. Coupling between the modes (energy transfer) in the horizontal 
direction occurs due to changes in sound speed, density, or bathymetry. The Chiu-Ehret 
model [Ref 32] determines the pressure field using numerical evaluations. Sagos [Ref 33] 
modified this model for three dimensions. 

Chiu [Ref 35] then developed a model in MATLAB™ called BBCM. In this model the 
spectral pressure at the receiver is decomposed by an FFT algorithm. At each point in the 
pressure field and at each frequency, the amplitude envelope for each existing normal mode 
is calculated, then coupled to determine the total pressure by 



V(r,zj) = f S(/)^-^U n (r;/)Z n (z;r,/)ex p{/(j Q k n {r'\f)dr- 2^/) 



■di (10) 



where S(f) is the source amplitude spectrum, Zn is a linear combination of the normal 
modes, and U t (r;j) is the amplitude slowly varying modulation complex envelope. It is in 
this calculation that the coupling of the normal modes occurs. 

In the MATLAB™ version, it has several advantages: 

• Can handle both discrete and broad band signals defined by S(f), 

• The bottom can be varied— not limited to a smooth slope, 

• FFT algorithms are extremely fast. Pressure fields are calculated in a fraction of 
the time that would be taken by URTEXT. 

BBCM uses a rigid bottom for the waveguide, while in the upcoming comparison no 
such bottom is assumed. Therefore, the effect of this rigid bottom on the sound in the water 
column must be examined. 
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C. THE TWO LAYER WAVE GUIDE 



Figure 3 shows a sound velocity profile of a hypothetical wave guide. In this guide the 
boundary conditions are a pressure release surface and a rigid bottom, with a discontinuous 
change in sound speed at depth L 




Figure 3 - Sound Velocity Profile for Comparison 

Normal mode theory produces an exact solution to the normalized Helmholtz equation 
for sound pressure P. This equation, given in cylindrical coordinates for a point source, is 



l_d_ 


d ' 


_ d 2 


(JO 




r — 


P+ — ^P + 




r dr 


dr. 


dz z 


c(z) 



P = -6(r)6(z-^)e 



JOit 



( 11 ) 



and has a solution of 



P(r,z,t) = -jne^'^Z^ZojZSz^Hf'iKj)] (12) 

n 

where Z n (z ) is the depth dependent normal mode, Zq is the depth of the source, r is the 
range from the source, K n is the horizontal propagation constant of the nlh normal mode. 
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and /^, 2) is the Hankel function where = J 0 - jY 0 [Ref 34, p. 431J. The function Z n (z) 
is a solution to the homogeneous Helmholtz equation 



d 2 , 

2 k.„Z = 0 

dz 2 " 2n " 



(13) 



First we examine the case where H- Z n (z) must satisfy the boundary conditions of 
Z(0)=0, continuity of particle velocity and pressure across the interface at z = L, and Z n {z) 
must approach zero as z goes to infinity. 

For the above conditions, we guess at two functions Z, n and Z^ n to represent Z„(z), 



An(z) = A n sin(k lzn z) 



0 <z<L 



(14) 



and 



Z ln (z) = B n e- fi - (z - L) L<z (15) 

where k lzn is the vertical wave number of the nth depth dependent normal mode in the upper 
layer and is the vertical propagation number of the exponential function in the lower 
layer. The numbers k lzn and P„ are related by the equation 





\ 2 


/ \ 






03 


= 


k l/iz ~ 


— 






\ c u 






(16) 



Evaluation of (14) at the surface proves that it is a valid solution, but the behavior of the 
functions at the interface of z = L needs to be determined. At the interface, continuity of 
sound pressure and normal particle velocity must hold. Therefore, 



A n S ' n { k Kn L )- B n 



(17) 



and 



( 18 ) 

Combining (17) and (18) and manipulating gives the transcendental equation 



12 



tan(*n^)=-7f J 

V>n L 



( 19 ) 



which, by defining 



y=Kz L 



and 



2 i 2 

a = L 



\ 2 
03 

\ C 1/ 



\ 2 

0) 

\ C l/ 



yields 



( 20 ) 



( 21 ) 



tan y = - 




( 22 ) 



This equation can be used to solve for k Un , and this yields p„and k„, when y < a. 
Normalization of Z n {z) then produces the terms A„ and B n , giving the final solution to (13). 

This solution corresponds to an oscillatory function above L overlying a function which 
decays exponentially as z increases. An example, Z^z), at a frequency somewhat above 
cutoff, is shown in Figure 4. As can be seen, the sound energy is trapped in the upper layer 
with some energy leaking into the bottom. 
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Figure 4 - Fluid Layer Overlying an Infinite Bottom 

The cutoff frequency of each normal mode can be determined at each point where a n 
approaches an asymptote of tan y in (22), or where 



x 3ji (2m - 1)ji 



a n 9 ’ o ’■ 



n - 1,2,3, ... 



(23) 



Substitution of (21) into (23) and manipulation yield a cutoff angular frequency of 
normal mode //, 






(2// - l)flr,c 2 
2Lj(c*-c, 2 ) 



(24) 
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We now examine the case where there is a rigid bottom at z-H. The solutions of the 
Helmholtz equation now must fulfill the additional boundary condition < 7—~ ■= 0 



dz 



Z-H 



For the revised boundary' conditions, we choose 



Z,„(z) = A,sin(/: Un z) 



0 < z < L 



(25) 



and 



Zj„(z) = fl„cosh[p n (//-z)] L < z < H (26) 

The hyperbolic cosine was chosen in analogy with the exponential choice in the previous 
case. It satisfies the boundary conditions at z = H. 

Next is evaluation of (25) and (26) at the interface of z = L. At the interface, 

aMK.L)- B n cosh[p„ (// - Z.)] (27) 

and 



AK, c °s (KJ-) - -BA sinh[p„(H - i)] 
Combining (27) and (28) gives 



(28) 



(29) 



which, by substitution of y and a yields 



tan y =- 



ycoth[(H/ L -i)J7~f 

4^7 



(30) 



This equation can be used to solve for k izn and then p„and K n can be obtained. When y<a 
the solution is similar to that for the infinite bottom condition in that it corresponds to an 
oscillatory function above Land decays as z increases from L to H. The sound energy in 
this case is trapped in the upper layer with some leakage into the lower layer (Figure 5). 
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Note that, from (30), as //decreases to L,y decreases to (n - \/2)n for n = 1, 2, 3,. . . 
and the normal modes will tend to approach those for a singe layer waveguide with a rigid 
bottom and a pressure release surface. 




Figure 5 - Two Layers Overlying a Rigid Sub-bottom 

The propagation of the normal mode exists in three regimes dependent upon k lzn . We 
determine (3,, in terms of k Un from (16) to be 

7 \ 2 / \ 2 ' 




U) 



\ c i / 



u> 

\ C 2j 




(31) 
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If c < 



(co) 


2 


(co) 


— 


— 


— 


\ c lJ 







, this is equivalent to y 
y > a then p„ becomes imaginary' and Zj n (z) 



< a, and therefore (25) through (30) hold. If 
becomes 



= ~B n co{k 2ln (H - z)] 
where k lzn is defined by 



(32) 



k 2 = lc 2 - 

K lzn K lzn 



,.y 



U) 



\ c l/ 



\ 2 



CO 



\ c lJ 



(33) 



and is related to K n by 

** = (-) "4, 



(34) 



By applying the boundary contains to (25) and (32) and making the appropriate 
substitutions for y and a we produce 



tan y 



y cot [(%~ 1 )V > ;2 ~ <7 






(35) 



This equation describes a normal mode trapped in the waveguide with a pressure release 
surface and a rigid bottom. How that mode propagates is dependent upon K n . 



co 



From (16), if — 

W 



< k* zn < I — j , then k n is real, and the mode propagates 



horizontally in the overall waveguide. If — < k* zn , K n becomes imaginary. The normal 



mode thus decays as it propagates horizontally. Figure 6 shows the vertical component of 
the normal mode which satisfies (35). 
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Normalized Amplitude Z 



Figure 6 - Two Layers Overlying a Rigid Sub-bottom 

The condition where p„ becomes imaginary is the major difference between the infinite 
bottom and the rigid, finite bottom. In the infinite bottom case, decreasing the layer depth L 
decreases the number of modes which can be propagated. However, the finite bottom adds a 
propagation mode throughout waveguide which must be corrected if an accurate comparison 
is to be made. 

Another difficulty encountered is that the above discussion is valid only when all 
interfaces arc parallel. With parallel layers, the grazing angle 0 is easily computed. With a 
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sloping interface L=f(r) this angle changes with range. It can be approximated numerically, 
but no analytic solution exists. 

To extend to a sloping interface in a waveguide consisting of two layers with a rigid 
bottom of depth H and fluid interface at depth L, all the normal modes with cutoff 
frequencies below the frequency of the source are excited and propagate according to (25) 
and (26). As the sound goes upslope, each normal mode approaches its first transition depth 
where (3 n <= 0 and the pressure amplitude of the exponential function becomes uniform over 

the entire lower layer (Figure 7). After the transition, the normal mode now behaves 
sinusoidally with depth in both layers (Figure 6). As L continues to decrease, the mode 
approaches the second transition point, where K n becomes imaginary and the mode 

becomes evanescent in the horizontal plane. 
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Figure 7 - Fluid Layer over a Rigid Bottom at First Transition 
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III. PROCEDURES AND RESULTS 



A. PROCEDURES 

For this research, the image model (URTEXT) and the coupled mode model (BBCM) 
were compared for the upslope case. Both models were run using MATLAB™, with 
URTEXT being translated from FORTRAN. Validation of this translation appears in 
Appendix A. 

The parameters chosen were the benchmark wedge problem of the Journal of the 
Acoustic Society of America [Ref 36] (Figure 8). The speed of sound in the water is 
uniform at 1500 m/s and the standard seawater density of 1.021 g/cm 3 is used. A point 
source transmitting a continuous wave at 25 Hz is located 4 km from the apex of the wedge 
and 100 m below the surface. The bottom slopes at 2.86°. The composition of the bottom 
will be varied to test different scenarios. 
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Three cases were run for comparisons. Scenario 1 uses a pressure release bottom. 
Scenario 2 has a fast bottom with no absorption loss, and Scenario 3 uses the same bottom 
with some absorption. Parameters for the environments are given in Table 1. 



TABLE 1 - ENVIRONMENTAL PARAMETERS 





Scenario 1 


Scenario 2 


Scenario 3 


ci (m/s) 


1500 


1500 


1500 


C2 (m/s) 


343 


1700 


1700 


Pi(g/cm 3 ) 


1.021 


1.021 


1.021 


p 2 (g/cm 3 ) 


0.00121 


1.500 


1.500 


a (dB/X) 


0.00 


0.00 


0.5 



One difficulty was the difference in measured values between the two models. BBCM 
uses SI units and produces transmission loss (among other values) in dB. URTEXT uses 
dimensionless parameters scaled to the dump distance determined by (8a) and produces a 
dimensionless pressure amplitude. This can be converted to dB by 

rL-201og^J (36) 

where TL is the transmission loss in dB, P is the dimensionless complex pressure, and Xc is 
the dimensionless dump distance. For each scenario, it was necessary' to convert the SI units 
into the dimensionless input units for URTEXT. Actual numbers will be given later. 

BBCM uses input files w'hich model the environment. These files define the overall 
waveguide (water column and underlying bottom) by describing the sound velocity profile, 
density profile, receiver ranges, receiver depths, attenuation of each normal mode, and source 
amplitude spectrum [Ref 35]. Depending upon the degree of resolution desired, these input 
files can get quite large. A program, SCENEGEN, was written to develop these input files 
specifically for the two dimensional wedge. 
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Another problem of BBCM concern the rigid bottom of the overall waveguide. As seen 
from the development in Chapter II this rigid waveguide floor affects the sound in the water. 
Depending on the overall depth of the waveguide, free propagating normal modes of the 
entire waveguide are excited which are coupled with the modes which are trapped only in the 
water layer. These free normal modes must be identified and eliminated for accurate results. 

This could be done by applying a range dependent attenuation coefficient for each 
normal mode, so the mode is attenuated by the bottom after it reaches its first transition 
point. However, BBCM does not support this. The attenuation vector that is entered in the 
model is range independent. Therefore, the effects of the rigid sub-bottom cannot be 
eliminated for this problem. 

B. SCENARIO 1 -PRESSURE RELEASE BOTTOM 

A pressure release bottom was simulated for both models. Since both models are 
incapable of running if either q or p 2 - 0, an air-water interface was deemed sufficient for 
the bottom. Inputs for URTEXT are given in Table 2. As the sound goes closer to the apex, 
adiabatic theory predicts less pressure as cutoff depth for each successive normal mode is 
reached. The actual outputs are seen in Figure 9. 



TABLE 2 - URTEXT INPUTS FOR COMPARISONS 





Scenario 1 


Scenario 2 


Scenario 3 


p 


3.0° 


2.86° 


2.86° 


Y 


1.43° 


1.43° 


1.43° 


y 

a 2 


4.412 


0.8824 


0.8824 


y 

/ P: 


850 


0.6667 


0.6667 




69.87 m 


638 m 


638 m 


r 1 


57.33 X c 


6.27 X c 


6.27 X c 


1 % 


0.0001 dB 


0.0001 dB 


0.07963 dB 
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For the inputs from Table 2, normal mode theory for a rigid bottom predicts seven 
modes will be excited by the source. But since the source is located at the pressure node of 
the even numbered modes, these modes are not excited. Therefor, only four modes are 
excited and cutoff points for these modes should be observed. The URTEXT results show 
these predicted characteristics, with cutoff points for the first mode at 3400 m, the third 
mode at 2200 m, and the fifth mode at 1000 m. Cutoff for the seventh mode cannot be seen 
due to the high density of contour lines in Figure 9. Further information is seen in the 
transmission loss plots of Figures 10 through 13. 

BBCM did not produce usable data for this scenario. This is not surprising, when one 
considers the original purpose of the model. BBCM was written for acoustic tomography 
calculations. As such, the pressure release case is not considered, since there are no actual 
pressure release bottoms in the ocean. The model ends up taking the sound from the water 
column and trapping it in the underlying bottom. 




0 260 520 780 1040 1300 1560 1820 2080 2340 2600 2860 3120 3380 3640 3900 



Range from Source (m) 



Figure 9 - Scenario 1 - URTEXT Transmission Loss in dB re 1 ni 
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Figure 10 - Scenario 1 - 30 m Receiver 




Figure 11 - Scenario 1 - 50 m Receiver Depth 
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Figure 12 - Scenario 1 - 100 m Receiver Depth 
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C. SCENARIO TWO -FAST, LOSSLESS BOTTOM 

In this scenario a more realistic environment was used. Inputs were as seen in Table 1. 
This produced the inputs for URTEXT as given in Table 2. Figures 14 through 17 are 
transmission loss profiles of both URTEXT and BBCM within the water column. 

Only the water column is considered as URTEXT does not calculate pressures in the 
underlying bottom. However, areas of significant difference can be seen. These are 
primarily due to the effect of the rigid sub-bottom. Areas where BBCM show an increase in 
transmission loss are likely due to a normal mode going through transition, while areas 
where there is significant decrease in transmission loss are possibly due to modal 
interference. 




Figure 14 - Scenario 2 (Lossless Bottom) - 30 m Receiver Depth 
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Figure 15 - Scenario 2 (Lossless Bottom) - 50 m Receiver Depth 




Range from Souroe (m) 



URTEXT 

BBCM 



Figure 16 - Scenario 2 (Lossless Bottom) - 100 m Receiver Depth 
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Figure 17 - Scenario 2 (Lossless Bottom) - 150 m Receiver Depth 



D. SCENARIO 3 - LOSSY BOTTOM 

The environment for this run was the same as for Scenario 2, except a loss term of 0.5 
dB/X was added for the bottom. This converted to 0.07963 dB for the loss term in 
URTEXT. Figures 18 through 21 show the results of this run. 

For this scenario, another difficulty was encountered. While the benchmark calls for 
attenuation by the bottom for the source frequency, that attenuation must be applied to the 
normal modes of interest. Therefore, a range dependent attenuation coefficient must be 
determined. Since no such coefficient can be applied, the results developed are rough 
approximations at best near or around modal cutoff. 

BBCM runs were not accurate but are estimates. BBCM was run with three normal 
modes. Application of the attenuation term is a difficult exercise, first because it needs to be 
range dependent to more accurately eliminate the rigid sub-bottom, and second because the 
term must be distributed among the modes. Further research in this problem is necessary. 
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Again, good correlation is noted between the models. The addition of the attenuation 
term has no great effect on the sound in the water column. However, there is a significant 
drop-off at the end of the BBCM graphs. These are the points where the receivers are now 
in the bottom. Here the attenuation is significant. 

Another significant feature is in Figure 18 at 150 m range. Here a dip is seen in 
URTEXT which is not seen in BBCM. At this relatively short range surface interference is 
noted. URTEXT picks this up (see Appendix A) but BBCM does not. Again, the original 
purpose of BBCM is to make predictions at relatively long ranges, significantly greater than 
the influence of surface interference [Ref 35]. At these longer ranges, BBCM is a robust 
model. 
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Figure 19 - Scenario 3 (Lossy Bottom) - 50 m Receiver Depth 
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Figure 20 - Scenario 3 (Lossy Bottom) - 100 m Receiver Depth 
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Figure 21 - Scenario 3 (Lossy Bottom) - 150 m Receiver Depth 
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IV. CONCLUSIONS AND RECOMMENDATIONS 



Figures 22 through 25 show the results of Jensen and Ferla using a PE method 
compared to the results of URTEXT [Ref 36]. 
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Figure 22 - Scenario 2 - 30 m comparison of PE and URTEXT 
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Figure 23 - Scenario 2 - 150 m Comparison of PE and URTEXT 
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Figure 25 - Scenario 3 - 150 m Comparison of PE and URTEXT 

For the cases examined— upslope propagation over both a lossless and lossy bottom, 
both URTEXT and BBCM have good agreement between the two. In the close range, 
BBCM does not show surface interference while URTEXT does. However, BBCM does 
develop a pressure field in the underlying bottom while URTEXT is effective only in the 
water column. Furthermore, BBCM is much quicker to create the pressure field. Using a 
range resolution of 25 m, BBCM created a pressure field for each scenario in 2.5 hours per 
run. Using the same equipment, URTEXT took 10 hours each run to develop the same 
resolution without calculating transmission loss in the bottom. BBCM therefore is more 
time effective. 

Other strengths of BBCM are its flexibility. Because the input files model of the 
environment, different conditions can be simulated by changing the input files. Upslope, 
downslope, cross-canyon, stratified water or bottom layers can all be simulated. URTEXT is 
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limited only to the condition of the wedge shaped bottom. It is still a good tool for analysis 
of the specific problem, but BBCM is a more flexible tool for real world conditions. 

Areas of further research present themselves. URTEXT needs to be sped up. 

MATLAB is very fast with matrix arithmetic and manipulation. A revision of URTEXT 
which takes advantage of this capability would be handy in the analysis of the wedge 
problem. Also, new developments to the image problem include using doublets and sets of 
four images. Perhaps including these as subroutines would speed up URTEXT. 

The examination of a sloping interface over a rigid sub-bottom needs to be investigated. 
A modification of URTEXT to include this condition promises interesting results. 
Furthermore, a solution of the above problem which take incident angle needs to be 
determined. 

Finally, both the downslope and cross-slope cases need to be investigated with BBCM. 
The upslope case is a standard comparison with other models, but predictions in either 
conditions might prove useful. 
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APPENDIX A 



To develop a comparison standard between the image model and the coupled mode 
model, we translated the image program URTEXT from FORTRAN to MATLAB™ [Ref 
37]. MATLAB is a high level scripting language which runs on a variety of computer 
operating systems including VAX, MS-DOS, UNIX, and Macintosh. This improved the 
flexibility of the image model to let it be run on desktop computers and advanced 
workstations. Furthermore, since the coupled mode model is also written in MATLAB, 
comparison of CPU run times can be made. 

Validation of the MATLAB version of URTEXT, called URTEXTBAT, was conducted 
by comparing the results of URTEXTBAT to those of URTEXT as run on the IBM 3033 
located at the Naval Postgraduate School. Table A-l shows the input conditions of the 
programs. Pressure amplitudes were computed for scaling distances from three to ten dump 
distances from the wedge apex and at receiver angles from the surface to the bottom. 



TABLE A-l 



p 


10° 


Y 


5° 


r l 


2 


Pl/P2 


2.0 


C]/C2 


0.5 


a/k2 


0.001 





0.0 



Tables A-2, A-3, and A-4 show the results of the comparison. The major error at the 
surface is due to the difference in precision of the different languages. FORTRAN uses 
single precision variables, accurate to six decimal places. MATLAB is precise to 13 places. 
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therefore at the small pressures near the surface the round-off errors are extreme. Excluding 
the errors at the boundaries, the average absolute difference of the field was 0 . 00206 . 

TABLE A-2 



FORTRAN T2 



6 


3 


4 


5 


6 


7 


8 


9 


10 


10 


0.00324 


0.00813 


0.0026 


0.00158 


0.00118 


0.00118 


0.00153 


0.0031 


9 


1.1393 


0.75924 


0.63381 


0.4346 


0.30078 


0.2086 


0.13697 


0.0337 


8 


2.13245 


1.17088 


1.16923 


0.80956 


0.56091 


0.38328 


0.20159 


0.14456 


7 


2.86145 


1.4393 


1.52998 


1.07776 


0.74823 


0.47362 


0.34154 


0.34534 


6 


3.24737 


1.6604 


1.66464 


1.20835 


0.84126 


0.5739 


0.56706 


0.5354 


5 


3.26776 


1.93661 


1.56257 


1.19399 


0.8313 


0.79333 


0.74299 


0.67374 


4 


2.97124 


2.06111 


1.25278 


1.05485 


0.91045 


0.9172 


0.82486 


0.72998 


3 


2.51693 


1.97233 


0.80655 


0.84476 


0.98048 


0.90308 


0.79069 


0.68846 


2 


2.39419 


1.65642 


0.51872 


0.72407 


0.84396 


0.74538 


0.64007 


0.55031 


1 


1.65191 


1.13991 


0.3643 


0.54176 


0.54487 


0.46663 


0.39299 


0.33296 


0 


0.75925 


0.48278 


0.14244 


0.15337 


0.14283 


0.1 1185 


0.08627 


0.06728 



TABLE A-3 

MATLAB™ Ij 



6 


3 


4 


5 


6 


7 


8 


9 


10 


10 


0.0003 


0.0008 


0.0003 


0.0002 


0.0001 


0.0001 


0.0002 


0.0003 


9 


1.1378 


0.7536 


0.6319 


0.4337 


0.3004 


0.2087 


0.1377 


0.0358 


8 


2.1323 


1.1668 


1.1681 


0.8093 


0.5612 


0.3845 


0.2038 


0.1422 


7 


2.8627 


1.4368 


1.5293 


1.0781 


0.7494 


0.4822 


0.3388 


0.343 


6 


3.25 


1.6586 


1.6652 


1.2094 


0.8436 


0.571 


0.5643 


0.5336 


5 


3.2716 


1.9359 


1.5637 


1.1955 


0.8286 


0.7897 


0.7412 


0.6726 


4 


2.9753 


2.0614 


1.2543 


1.0564 


0.9016 


0.9152 


0.824 


0.7296 


3 


2.5161 


1.9736 


0.8081 


0.8451 


0.9778 


0.9026 


0.7908 


0.6887 


2 


2.3904 


1.6583 


0.5202 


0.721 


0.8439 


0.746 


0.6408 


0.551 


1 


1.6515 


1.1421 


0.3659 


0.5405 


0.5463 


0.468 


0.3941 


0.3339 


0 


0.7354 


0.485 


0.1428 


0.1557 


0.1449 


0.1134 


0.0875 


0.0683 
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TABLE A-4 



Absolute Difference t"2 


6 


3 


4 


5 


6 


7 


8 


9 


10 


10 


0.00294 


0.00733 


0.00230 


0.00138 


0.00108 


0.00108 


0.00133 


0.00280 


9 


0.001 50 


0.00564 


0.00191 


0.00090 


0.00038 


0.00010 


0.00073 


0.00210 


8 


0.0001 5 


0.00408 


0.00113 


0.00026 


0.00029 


0.00122 


0.00221 


0.00236 


7 


0.00125 


0.00250 


0.00068 


0.00034 


0.001 1 7 


0.00858 


0.00274 


0.00234 


6 


0.00263 


0.00180 


0.00056 


0.00105 


0.00234 


0.00290 


0.00276 


0.00180 


5 


0.00384 


0.00071 


0.00113 


0.00151 


0.00270 


0.00363 


0.00179 


0.00114 


4 


0.00406 


0.00029 


0.00152 


0.00155 


0.00885 


0.00200 


0.00086 


0.00038 


3 


0.00083 


0.00127 


0.00155 


0.00034 


0.00268 


0.00048 


0.0001 1 


0.00024 


2 


0.00379 


0.00188 


0.00148 


0.00307 


0.00006 


0.00062 


0.00073 


0.00069 


1 


0.00041 


0.00219 


0.00160 


0.00126 


0.00143 


0.00137 


0.00111 


0.00094 


0 


0.02385 


0.00222 


0.00036 


0.00233 


0.00207 


0.00155 


0.00123 


0.00102 



To see if the MATLAB is more precise, a run was made on an 18° slope with the data 
taken from the surface to 2° below sampled every 0.1°. Figure A-l shows pressure falls off 
linearly as the sound nears the surface. This implies the model is valid at the boundaries. 
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Receiver Range 



4 — 2 3 — ♦ — 4 □ 5 




Figure A-l - Boundary Condition Check 

To further explore the differences between the two models, a third model, called 
WEDGESLO was used. This model was validated by Kaswandi [Ref 14]. This program is 
written in BASIC to be run on an MS-DOS based computer, but used double precision 
variables. A single run at r?=5 showed correlation between WEDGESLO and 
URTEXTBAT to four places, including at the boundary. 

Some special cases were run on URTEXTBAT before it was considered valid. First, to 
verify reciprocity of the field, the source was set at r;=10 and y varied from 0° to 10°. The 
receiver was fixed at 6 = 5° and at r?= 2. Table A-5 shows the results and an error 
comparison. The errors are within the round off tolerance of the program, showing 
URTEXT as valid for both the upslope and downslope fields as well as v alidating the 
acoustic reciprocity of the model. 
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TABLE A-5 - RECIPROCITY CHECK 



6 


Pressure 


Difference 


10 


0 


0.00310 


9 


0.03349 


0.00021 


8 


0.14344 


0.00112 


7 


0.34385 


0.00149 


6 


0.53402 


0.00138 


5 


0.67274 


0.00100 


4 


0.72949 


0.00049 


3 


0.6885 


0.00004 


2 


0.55086 


0.00055 


1 


0.33387 


0.00091 


0 


0.06837 


0.00109 



Another test was run using matched acoustic impedance of the water and the bottom. 
The bottom then becomes transparent and the pressure field will be that of a surface 
interference effect (Lloyd's Mirror). Density and speed of sound ratios were set to .9999 
and slope angle was set to 18°. The receiver was set at the same depth as the source and 
moved downslope. The pressure field predicted is approximated by the equation 

pM)- — 

r 

where h and d are depths of the receiver and the source, A is the pressure amplitude at the 
source at one meter, and k is the wave number. Frequency used was 1000 Hz. 

Figure A-2 shows a plot of pressure as a function of range from the source. 
Interference maxima and minima from URTEXTBAT match those predicted by equation 



khd 



sin- 



(A-l) 



(A-l). 
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Figure A-2a - Close Range Surface Interferance 




Lloyd's Mirror - - - URTEXTBAT 



Figure A-2b - Long Range Surface Interferance 
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As seen from the above, URTEXTBAT has undergone extensive proofing. When 
compared with previous versions in FORTRAN and BASIC, it is more flexible both in 
platforms run and model variables. While the program does not take full advantage of the 
features of MATLAB, and has significant room for optimizations, the program in its present 
form is more than sufficient for the purposes of this research. 
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APPENDIX B 



%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% URTEXT for MATLAB 
% 

% Adapted from FORTRAN 
% 

%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Enter Input Variables 

% Required are: Range of Source and Receiver 
% Wedge angle in Degrees 

% Source and Reciever angles from the Bottom 
% Speed of sound ratio and density ratios 
% Bottom loss coefficient and cross slope range 
%%%%%%%%%%%%%%%%%%%%%%%% 

B=input( 'Enter bottom wedge angle: '); 

G=input( 'Enter source angle from the surface: '); 

D=input( 'Enter receiver angle from the surface: ' ); 

rhol=input( 'Enter Density ratio, water to bottom: '); 

CC=input( 'Enter speed ratio, water to bottom: '); 

rl=input( 'Enter range of source from apex: '); 

r2=input( 'Enter range of receiver from apex: '); 
AL=input( 'Enter bottom loss coefficient: '); 

yO=input( 'Enter cross-slope range: '); 

% Determine number of image pairs 

% and convert to radians. Also determine constants 

% Function rads (A) converts degrees into radians 
% Contained as a seperate M-File in the MATLAB Toolbox 

function X=rads(A) 

% Converts degrees into radians 

X=A*pi/180 ; 

end 



Nl=f ix ( 1 80/B ) ; 

B=rads ( B ) ; 

G=rads(G) ; 

D=rads(D) ; 
c2=CC^2 ; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Determine scaling for fast or slow bottom 
%%%%%%%%%%%%%%%%%%%%%%%% 

tqq=tan ( B ) ; 
if CC<1, 

tqql=acos ( CC ) ; 
tqq2=sin( tqql ) ; 
else 

tqql=acos ( 1/CC) ; 
tqq2=tan ( tqql ) ; 
end; 
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tqq3=2 * tqq2 * tqq ; 
t4=pi/tqq3 ; 
ql=l/ sqrt ( 2 ) ; 
d2=rl ^2+r2 ^2+y0~2 ; 
r3=2 *rl *r2 ; 
sl = l ; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Determine range to receiver for each image 
%%%%%%%%%%%%%%%%%%%%%%%% 

for n=l : 1 :Nl , 
if sl>0 

tl ( n ) = ( n-1 ) *B+G; 
else 

tl ( n) =n*B-G; 
end; 
s l=-sl ; 

r8(n)=sqrt(d2-r3*cos(tl(n)-D) ) ; 
r9(n)=sqrt(d2-r3*cos( tl(n)+D) ) ; 
end; 

pl=0;p2=0; 
for n=l:l:Nl, 

s2=(-l)^(fix(n/2) ) ; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Calculate Reflection Coefficient for Upper Space 
%%%%%%%%%%%%%%%%%%%%%%%% 

wl=2*c2*AL; 
il=fix( (n-1 ) /2 ) ; 
for i=l:l:il, 

s ( i ) =abs ( rl *sin( tl ( n) -2 *i*B ) +r2 * sin( 2 *i*B-D ) ) /r8 ( n) 
if s(i)>=l, 
s(i)=l; 
end; 

c(i)=sqrt(l-(s(i) A 2) ) ; 
t=s(i) /rhol; 
w0=( -c2+ ( c ( i ) ^2 ) ) ; 
y=sqrt( ( w0~2 ) + ( wl ^2 ) ) ; 
z=abs( wO) ; 
if y<=z, 
y=z; 
end; 

yl=ql*sqrt(y+wO) ; 
y2=-ql*sqrt(y-w0) ; 
zl=t-y2 ; 
z2=-yl ; 

z3=zl/ ( zl A 2+z2~2) ; 
z4=-z2 / ( zl"2+z2^2 ) ; 
zl=t+y2; 
z2=yl ; 

z5=zl*z3-z2*z4 ; 
z6=zl*z4+z2*z3; 
e(i)=z5; 
f (i)=z6; 
end; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Determine Pressure Contribution for each upper image 
%%%%%%%%%%%%%%%%%%%%%%%% 
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zl=0 ; z2=0 ; z3=0 ; z4=0; z5=l ; z6=0 ; 
if n>2, 

for i=l:l:il, 
zl=e ( i ) ; 
z2=f ( i ) ; 
z3=z5; 
z4=z6 ; 

z5=zl*z3-z2*z4; 
z6=zl*z4+z2*z3; 
end; 
end; 
zl=z5 ; 
z2=z6; 
t=t4*r8(n) ; 
z3=cos( t) ; 
z4=-sin( t ) ; 
z5=zl*z3-z2*z4 ; 
z6=zl*z4+z2*z3; 
pl=pl+s2*z5/r8 (n) ; 
p2=p2+s2*z6/r8 (n) ; 
il=il+l ; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Determine Reflection Coefficient for the Lower Images 
%%%%%%%%%%%%%%%%%%%%%%%% 

for i=l : 1 : il , 

s(i)=abs(rl*sin( tl(n)-2* (i-1) *B ) +r2*sin( 2* (i-1) *B+D) )/r9(n) 
if s(i)>l / 
s(i)=l; 
end; 

c ( i ) =sqrt ( 1.0001-s(i)^2) ; 
t=s( i) /rhol ; 
w0=-c2+c ( i ) ^ 2 ; 
y=sqrt ( wO ^2+wl ^2 ) ; 
z=abs(wO) ; 
if y<=z, 
y=z; 
end; 

yl=ql*sqrt(y+wO) ; 
y2=-ql * sqrt ( y-wO ) ; 
z l=t-y2 ; 
z2=-yl; 

z3 = zl/( zl /s 2 + z2 /v 2) ; 
z4=-z2/(zl"2+z2"2) ; 
zl=t+y2 ; 
z2=yl ; 

z5=zl *z3-z2*z4 ; 
z6=zl*z4+z2*z3; 
e(i)=z5; 
f (i)=z6; 
end; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Find Pressure Contribution for Lower Images 
%%%%%%%%%%%%%%%%%%%%%%%% 

zl=0 ; z2=0 ; z3=0 ; z4=0 ; z5=l ; z6=0 ; 
for i=l:l:il, 
zl=e( i ) ; 
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z2=f (i) ; 
z3=z5 ; 
z4=z6 ; 

z5=zl*z3-z2*z4; 
z6=zl*z4+z2*z3; 
end; 
zl=z5 ; 
z2=z6; 
t=t4*r9(n) ; 
z3=cos( t) ; 
z4=-sin( t) ; 
z5=zl*z3-z2*z4; 
z6=zl*z4+z2*z3; 
pl=pl+s2*z5/r9(n) ; 
p2=p2+s2*z6/r9(n) ; 
end; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Total Pressure at the reciever is; 
%%%%%%%%%%%%%%%%%%%%%%%% 

pz=sqrt (pl^2+p2' s 2 ) ; 
end; 
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%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% URTEXTBAT for MAT LAB 

% 

% Adapted from FORTRAN 

% 

%%%%%%%%%%%%%%%%%%%%%%%% 

% This program determines a 2 dimensional 
% Sound pressure field over a sloping bottom 
%%%%%%%%%%%%%%%%%%%%%%%% 

% Enter Input Variables 
% Required are: Range of Source 
% Maximum range of the Receiver 
% Wedge angle in Degrees 
% Source angle from the Surface 
% Speed of sound ratio and density ratios 
% Bottom loss coefficient and cross slope range. 

% The program determines receiver angles over 10 increments 
% and range from source range + one dump distance to the 
% maximum receiver range 
%%%%%%%%%%%%%%%%%%%%%%%% 

B=input( ’Enter bottom wedge angle: '); 

G=input( 'Enter source angle from the surface: '); 

rhol=input( 'Enter Density ratio, water to bottom: '); 

CC=input( 'Enter speed ratio, water to bottom: '); 

rl=input( 'Enter range of source from apex: '); 

range=input( 'Enter maximum receiver range: '); 

AL=input( 'Enter bottom loss coefficient: '); 

b=B / 1 0 ; countk=0 ; countl=0 ; 

y0=input( 'Enter cross slope range: '); 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Determine number of image pairs 
% and convert angles to radians 
%%%%%%%%%%%%%%%%%%%%%%%% 

time=clock; 

Nl=fix( 180/B) ; 

B=rads(B) ; 

G=rads(G) ; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Loop for field calculation 
%%%%%%%%%%%%%%%%%%%%%%%% 

for k=rl + l :range/10 -.range, 
r2=k; 

countk=countk+l ; 
for 1=0 :b :b* 10 , 

D=1 ; countl=countl+l ; 

D=rads ( D) ; 
c2=CC / '2; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Find scaling distance for fast or slow bottom 
%%%%%%%%%%%%%%%%%%%%%%%% 

tqq=tan( B ) ; 
if CC<1, 
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tqql=acos(CC) ; 
tqq2=sin( tqql ) ; 
else 

tqql=acos( 1/CC) ; 
tqq2=tan ( tqql ) ; 
end; 

tqq3=2 *tqq2 *tqq; 
t4=pi/tqq3; 
ql=l/ sqrt ( 2 ) ; 
d2=r 1 ^2+r2 ^2+y(K2 ; 
r3=2 *rl*r2 ; 
sl=l ; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Determine range for each image to reciever 
%%%%%%%%%%%%%%%%%%%%%%%% 

for n=l:l:Nl / 
if sl>0 

tl (n) =( n-1 ) *B+G; 
else 

tl(n)=n*B-G; 
end; 
s 1 = — s 1 ; 

r8(n)=sqrt(d2-r3*cos(tl(n)-D) ) ; 
r9(n)=sqrt(d2-r3*cos( tl(n)+D) ) ; 
end; 

pl=0 ; p2=0 ; 
for n=l:l:Nl, 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Find Reflection coefficient for upper images 
%%%%%%%%%%%%%%%%%%%%%%%% 

s2=( -1 ) ^ ( f ix( n/2 ) ) ; 
wl=2*c2*AL; 
il=fix( (n-1 ) /2 ) ; 
for i=l:l:il, 

s ( i ) =abs ( rl *sin( tl ( n) -2 *i*B ) +r2 *sin( 2 *i*B-D ) ) /r8 ( n) ; 
if s(i)>=l / 

s(i)=l; 

end; 

c ( i ) =sqrt (l-(s(i)' s 2) ) ; 
t=s( i) /rhol; 
w0=(-c2+(c(i)^2) ) ; 
y=sqrt( ( w(K2 ) + ( wl ^2 ) ) ; 
z=abs( wO ) ; 
if y<=z, 

Y=z; 

end; 

yl=ql*sqrt (y+wO ) ; 
y2=-ql*sqrt(y-w0) ; 
zl=t-y2 ; 
z2=-yl ; 

z3=zl/(zl^2+z2^2); 
z4=-z2 / ( zl^2+z2^2 ) ; 
zl=t+y2 ; 
z2=yl ; 

z5=zl *z3-z2*z4 ; 
z6=zl*z4+z2*z3; 
e( i)=z5; 
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f (i)=z6; 
end; 

zl=0; z2=0 ; z3=0 ; z4=0 ; z5=l ; z6=0 ; 
if n>2, 

for i=l:l:il, 
zl=e( i ) ; 
z2=f (i) ; 
z3=z5; 
z4=z6 ; 

z5=zl*z3-z2*z4; 

z6=zl*z4+z2*z3; 

end; 

end; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Find Pressure contribution of upper images 
%%%%%%%%%%%%%%%%%%%%%%%% 

zl=z5 ; 
z2=z6 ; 
t=t4*r8(n) ; 
z3=cos( t) ; 
z4=-sin( t ) ; 
z5=zl *z3-z2*z4 ; 
z6=zl*z4+z2*z3; 
pl=pl+s2*z5/r8(n) ; 
p2=p2+s2 *z6/r8 ( n) ; 
il=il+l ; 
for i=l:l:il, 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Find refection coefficient for lower images 
%%%%%%%%%%%%%%%%%%%%%%%% 

s ( i ) =abs ( rl *sin( tl ( n) -2 * ( i-1 ) *B ) +r2 *sin(2*(i-l) *B+D) ) /r9 ( n ) ; 
if s(i)>l, 
s ( i ) =1 ; 

end; 

c(i)=sqrt( 1.0001-s(i) A 2); 
t=s( i) /rhol ; 
w0=-c2+c(i) ^2; 
y=sqrt(w0~2+wl~2 ) ; 
z=abs ( wO ) ; 
if y<=z, 

y=z; 

end; 

yl=ql*sqrt(y+wO) ; 
y2=-ql *sqrt ( y-wO ) ; 
zl=t-y2 ; 
z2=-yl ; 

z3=zl/(zl^2+z2~2); 
z4=-z2/( zl~2+z2~2) ; 
zl=t+y2 ; 
z2=yl ; 

z5=zl*z3-z2*z4; 
z6=zl*z4+z2*z3; 
e(i)=z5; 
f (i)=z6; 
end; 

zl=0 ; z2=0 ; z3=0 ; z4=0; z5=l ; z6=0; 
for i=l:l:il. 



50 



zl=e ( i ) ; 
z2=f ( i ) ; 
z3=z5; 
z4=z6 ; 

z5=zl*z3-z2*z4; 
z6=zl *z4+z2*z3 ; 
end? 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Find Pressure contribution for lower images 
%%%%%%%%%%%%%%%%%%%%%%%% 

zl=z5; 
z2=z6 ; 
t=t4*r9 ( n) ? 
z3=cos ( t ) ; 
z4=-sin( t ) ; 
z5=zl*z3-z2*z4; 
z6=zl*z4+z2*z3? 
pl=pl+s2*z5/r9(n) ? 
p2=p2+s2*z6/r9(n) ; 
end; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Total pressure at the reciever is: 

%%%%%%%%%%%%%%%%%%%%%%%% 

pz ( count 1 , countk ) =sqrt ( pi ^ 2+p2 *2 ) ; 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Reset variables 
%%%%%%%%%%%%%%%%%%%%%%%% 

clear c;clear c2;clear d2;clear e? clear f;clear i; clear il; 
clear n;clear pl;clear p2?clear ql?clear r3;clear r8? clear r9 
clear s?clear sl?clear s2?clear t?clear tl? 

clear tqq;clear tqql?clear tqq2;clear tqq3;clear wO; clear wl ; 
clear y?clear yl;clear y2;clear z?clear zl;clear z2; clear z3? 
clear z4;clear z5;clear z6? 
end; 

countl=0; 

e t=e time ( clock , time ) 
end; 
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%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% SCENE GEN . M 
% 

%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% This program creates input files for BBCM, 

% the coupled mode model by Dr. C.S. Chiu. 

% This program creates the inputs for a wedge shaped 
% ocean overlying a fluid bottom. 

% Both water and bottom are considered uniform in characteristics 
% in depth. No irregularities in the bottom depth are 
% generated, nor stratifications in sound velocity profile. 

% 

%%%%%%%%%%%%%%%%%%%%%%%% 

% Input all the pertinent data 
B=input( 'Enter wedge angle: '); 

B=rads(B) ; 

z=input( 'Enter maximum depth of the water column: '); 

w=input( 'Enter length of the waveguide: '); 

dw=input( 'Enter range resolution: '); 

dz=input( 'Enter depth resolution: ’); 

cl=input( 'What is the speed of sound in the water? '); 

c2=input( 'What is the speed of sound in the bottom? '); 

rl=input( 'What is the density of the water? '); 

r2=input( 'What is the density of the bottom? '); 

n=input( 'How many modes? ' ) ; 

ak=input( 'What is the bottom loss coefficient?'); 
f=input( ’What is the source frequency? '); 

% Build data files — range 
r=0 : dw: w; 

save r.dat r /ascii; 

% receiver depth 
zr=0 : 1 0 : z; 

save zr.dat zr /ascii; 

% attenuation matrix 
for i=l : 1 : n, 

attenu(n)=ak; 

end; 

save attenu.dat attenu /ascii; 

% sound speed profile 
a=w/dw; b=z/dz; 
t=tan(B) ; 
for i=l : 1 : a+1 , 
for j=l:l:b+l, 

if t>( j-1 ) *dz/(w-(i-l ) *dw) , 
c( j,i)=cl; 
else 

c( j,i)=c2; 
end; 
end; 
end; 

dummy=c ( 1 : b , : ) ; 
for i=l:l:b, 

dummy ( i , 1 ) =c2 ; 
end; 
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for i=2:l:a, 

dummy ( : , i)=dummy( : ,i-l ) ; 
end; 

c= [ c ; dummy ; dummy ; dummy ; dummy ; dummy ; dummy ; dummy ; dummy ; dummy ] ; 
save c.dat c /ascii; 

% density profile 
for i=l:l:a+l, 
for j=l:l:b+l, 

if t>( j-1 ) *dz/(w-(i-l ) *dw) , 
rho ( j , i ) =rl ; 
else 

rho ( j / i ) =r2 ; 
end; 
end; 
end; 

dummy=rho( 1 :b, : ) ; 
for i=l:l:b, 

dummy ( i, 1 )=r2 ; 
end; 

for i=2:l:a, 

dummy ( : , i ) =dummy ( : , i-1 ) ; 
end; 

rho= [ rho ; dummy ; dummy ; dummy ; dummy ; dummy ; dummy ; dummy ; dummy ; dummy ] ; 
save rho.dat rho /ascii; 

% Source Spectrum Density Profile 
ssp=[0 0;f-(K005 l;f l;f+0.005 1;1000 0]; 
save ssp.dat ssp /ascii; 
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